library(boot)

setwd("code/")
source("insurance_source_score_1_prep.R")

plot_slopes_return_lm_for_boot <- function(
                                         d,
                                         w,
                                         propensity="MARKET",
                                         outcome="FAVOR>=3"
                                         ) {

    ## note - propensity is not the propensity score, it is training variable
    d$propensity <- d[,propensity]
    d$outcome <- d[,outcome]

    propensity_out <- lm(
        propensity ~ poly(EDUC, 2) + poly(INCOME, 3) +
            poly((AGEN), 5) +BLACK+HISP+ASIAN+MALE+RETIRED,
        data=subset(
            d[w,],
            DATN > 60 & !is.na(AGEN) & !is.na(INCOME) & !is.na(EDUC)
        ),
        weights = subset(
            d[w,],
            DATN > 60 & !is.na(AGEN) & !is.na(INCOME) & !is.na(EDUC)
        )$WEIGHT
    )


    d$PSCORE3 <- predict(
        propensity_out,
        newdata=d
    )

     lm_out_all <- lm(
        outcome ~ PSCORE3*I(PERIOD > 8.5),
        data=d[w,],
        weights = d[w,]$WEIGHT
    )

    return(
        summary(lm_out_all)$coef[which(rownames(summary(lm_out_all)$coef)=="PSCORE3:I(PERIOD > 8.5)TRUE"),1]
    )

}

sdat.under65.2$FAVOR01 <- sdat.under65.2$FAVOR >= 3
sdat.under65.2$COVERED0 <- sdat.under65.2$COVERED == 0


set.seed(987654321)                     #making this replicable

the_boot_uninsured <- boot(data = subset(
    sdat.under65.2,
    INCOME>40
    & !is.na(WEIGHT)
    ),
  statistic = plot_slopes_return_lm_for_boot,
  outcome="FAVOR01",
  propensity = "COVERED0",
  R = 10000,
  stype = "i"
  )

the_boot_market <- boot(data = subset(
    sdat.under65.2,
    INCOME>40
    & !is.na(WEIGHT)
    ),
  statistic = plot_slopes_return_lm_for_boot,
  outcome="FAVOR01",
  propensity = "MARKET",
  R = 10000,
  stype = "i"
  )


the_boot_selfinsure <- boot(data = subset(
    sdat.under65.2,
    INCOME>40
    & !is.na(WEIGHT)
    ),
  statistic = plot_slopes_return_lm_for_boot,
  outcome="FAVOR01",
  propensity = "SELFINSURE",
  R = 10000,
  stype = "i"
  )


print(the_boot_market$t0)
print(round(boot.ci(the_boot_market, type="perc")$percent[4:5], 2))
print(the_boot_uninsured$t0)
print(round(boot.ci(the_boot_uninsured, type="perc")$percent[4:5], 2))
print(the_boot_selfinsure$t0)
print(round(boot.ci(the_boot_selfinsure, type="perc")$percent[4:5], 2))
